Generalized scattering-matrix approach for magneto-optics in periodically patterned 

multilayer systems 
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We present here a generalization of the scattering-matrix approach for the description of the 
propagation of electromagnetic waves in nanostructured magneto-optical systems. Our formalism 
allows us to describe all the key magneto-optical effects in any configuration in periodically patterned 
multilayer structures. The method can also be applied to describe periodic multilayer systems 
comprising materials with any type of optical anisotropy. We illustrate the method with the analysis 
of a recent experiment in which the transverse magneto-optical Kerr effect was measured in a Fe film 
with a periodic array of subwavelength circular holes. We show, in agreement with the experiments, 
that the excitation of surface plasmon polaritons in this system leads to a resonant enhancement of 
the transverse magneto-optical Kerr efi^ect. 

PACS numbers: 78.20.Bh, 78.20.Ls, 78.66.Bz, 78.40.Kc 



I. INTRODUCTION 

In recent years a lot of attention has been paid to the 
study of the optical properties of nanostructured materi- 
als with both plasmonic and magneto-optic activity.^ The 
key idea is to use hybrid nanostructures containing both 
metals, which exhibit plasmon resonances, and ferromag- 
netic materials, which provide high magneto-optical ac- 
tivity for reasonably low values of the applied magnetic 
field, to profit from the best of the worlds of plasmon- 
ics and magneto-optics. In the context of these hybrid 
structures there are two main topics of interest. The first 
one is the use of the localization of the electromagnetic 
field due the excitation of the plasmons supported by the 
free electrons in metals to enhance the magneto-optical 
signals (Kerr effect, Faraday effect, etc.)<^i^ The nanos- 
tructuring in these magneto-plasmonic structures plays 
a fundamental role for several reasons. First of all, it 
provides a convenient way to couple the light of an ex- 
ternal source to the plasmons supported by these hybrid 
systems, avoiding so the typical wave vector mismatch in 
unstructured systems. On the other hand, by nanostruc- 
turing these hybrid systems one can manipulate light at 
the nanometer scale in several ways. In particular, the 
enhancement of the electromagnetic field can be largely 
increased since one can concentrate light in reduced vol- 
umes. This can be done either by making use of local- 
ized plasmon resonances in isolated structures such as 
wireSf^!^ diska^"— , or particles j^iiS or by a periodic perfo- 
ration in an otherwise continuous filmiiir— 

A second topic of interest is the use of magneto-optical 
effects to externally control either the properties of the 
transmission through perforated membranesi^r— or the 
very value of the surface plasmon wave vector j^i^^— In 
this latter case, the relevant configuration is the trans- 
verse magneto-optical Kerr effect (see Sec. IIIII below) 



since the other configurations induce a polarization con- 
version that implies a decoupling of the plasmon. 

In view of the relevance of these novel hybrid struc- 
tures, and in order to guide their design, it is cru- 
cial to have theoretical methods that are able to de- 
scribe the wave propagation in nanostructured magneto- 
plasmonic systems. A powerful approach, which is 
widely used to describe nanostructured systems with- 
out magneto-optical activity, is the so-called scattering- 
matrix formalismi^^^— In recent years, this approach has 
been extended to study different magneto-optical effects 
in nanostructured multilayer structures^* and to describe 
the wave propagation in periodic structures containing 
certain types of anisotropic mediai^ However, there are 
still basic physical situations which lie out of the scope 
of the existent implementations of the scattering formal- 
ism. Thus for instance, the Kerr and the Faraday effects 
in the transverse configuration, in which the magnetic 
field (or the magnetization of the sample) is parallel to 
the sample plane but perpendicular to the plane of inci- 
dence, cannot be addressed with the existent scattering- 
matrix-based approaches. More generally, when the op- 
tical anisotropy of the the materials involves off-diagonal 
elements of the permittivity tensor along the growth di- 
rection of the multilayer structure, none of the existing 
implementations of the scattering-matrix approach can 
describe the wave propagation in such structures. The 
technical problem lies in the fact that in such situations 
the propagating eigenstates in the different layers cannot 
be simply described with a standard eigenvalue problem. 
In this work, we show how this problem can be solved 
and present a generalization of the scattering-matrix ap- 
proach to describe the magneto-optics of hybrid nanos- 
tructured systems in any configuration. Moreover, the 
method can be applied to periodically patterned multi- 
layer systems comprising any kind of optically anisotropic 
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materials. We illustrate the capabilities of our formalism 
by addressing a recent experiment in which the transverse 
magneto-optical Kerr effect (TMOKE) was measured in a 
periodically perforated Fe film.^"* We show that, in agree- 
ment with the experimental results, the excitation of sur- 
face plasmon polaritons in these structures leads to an 
enhancement of the TMOKE signal. More importantly, 
our theoretical method paves the way for studying the 
interplay between plasmon-driven effects and magneto- 
optics in a wide variety of hybrid nanostructures. 

The rest of the paper is organized as follows. In sec- 
tion [n] we explain in detail how the scattering-matrix 
approach can be generalized to describe the wave prop- 
agation in all kind of magneto-optical and optically 
anisotropic periodic multilayer systems. Section lllll is de- 
voted to the analysis of the experiments of Ref. 14, which 
allows us to illustrate the power of the method. Then, 
we shall summarize the main conclusions of our work in 
section ITVl Finally, several technical issues related to the 
formalism developed in section HIl are addressed in detail 
in three appendixes. 



II. GENERALIZED SCATTERING-MATRIX 
APPROACH 

Our central goal is to solve the Maxwell's equations 
for a patterned multilayer structure containing any com- 
bination of materials (isotropic and anisotropic). For this 
purpose, we shall generalize the scattering-matrix ap- 
proach developed by Whittaker and Culshaw in Ref. |2^. 
Following this work, we shall first discuss the Maxwell's 
equations to be solved. Then, we shall address the band 
structure of an unbounded layer to determine the propa- 
gating eigenstates in the different layers. Then, we shall 
discuss how to construct the fields in a multilayer struc- 
ture using those eigenstates and, finally we shall describe 
how the scattering matrix can be used to determine the 
field amplitudes in the whole structure. 



A. Maxwell's equations 

Let us start by describing the Maxwell's equations 
to be solved. Assuming a harmonic time dependence 
exp(— iwt), the Maxwell's equations adopt the following 
form: V • eocE = 0, V-H = 0, VxH = -iweoeE, and 
V X E = iw/LtoH, where the permittivity is in general a 
tensor given by 



We now introduce the following rescaling: weqE — ^ E and 
^JJIqCquj = bj /c — > w. Thus, the final two equations to be 
solved are 



V X H = -ieE, 

V X E = iw^H. 



(2) 
(3) 



We consider here multilayer systems in which each 
layer can be, in principle, periodically structured. Thus, 
the tensor e is independent of z, where z corresponds to 
the growth direction of the structure, and it depends on 
the in-plane coordinates r = (x, y) in a periodic fash- 
ion. Due to this periodicity, it is convenient to work in a 
momentum representation for the in-plane coordinates. 
Thus, for a given Bloch wave vector k, we can expand 
the fields as a sum over reciprocal lattice vectors G 

H(r,z) = ^Hk(G,z)e^(''+«)-, (4) 

G 

E(r,z) = ^Ek(G,z)e*(''+°)•^ (5) 



Following Ref. |25|, we define the Fourier space vectors 

T 



h{z) = ^Hk(Gi,z),Hk(G2,z),.. 
e(z) = [Ek(Gi,z),Ek(G2,z),.. 



(6) 
(7) 



Note that, although Hk and Ek depend on k, the whole 
calculation is carried out for a fixed k, so such labels will 
be omitted in other symbols. 

In what follows, we shall need the Fourier expansion 
of the components of the permittivity tensor for the dif- 
ferent layers 



KG) 



unit cell 



dv ejj(r)e 



-iGr 



(8) 



where j, j — x, y, z, S is the area of the in-plane unit cell, 
and the matrix is such that (eyOcG' = — G'). 

Analogously, the components of the index tensor rjij (r) — 
[e~^{r)]ij have Fourier expansions fjij(G) and matrix rep- 
resentations fjij . 

With the notation just introduced, the momentum rep- 
resentation of a product such as e^ E becomes e^j-e. Thus, 
the relevant Maxwell's equations, Eqs. ([2]) and ([3]), can 
be written, in component form, as 



ikyh^z) - h'yiz) = -i'^e^.jejiz) 
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(9) 



Kiz)-iKK{z) = -ij^^yj^jiz) (10) 



(1) 



<^zy 



Notice that we have assumed that fl — 1 since we are 
interested in the optical regime. Notice also that for har- 
monic fields, the first Maxwell's equation is implied by 
the third one and the second can be satisfied by expand- 
ing the magnetic field in basis states with zero divergence. 



and 



ikxhy(^z) ikyhx{z^ 



ikyCziz) - e'y{z) 
e'xiz) - ik^Bziz) 
ikxGyi^z^ ikyex{z) 



lJ2^zjej{z), (11) 



= iiJ^hx{z) 
= iLo^hy{z) 
= iuj'^hziz), 



(12) 
(13) 
(14) 
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where kx and ky are diagonal matrices with {kx)GG = 
{kx + Gx) and (fcj,)GG = {ky+Gy), and the primes denote 
differentiation with respect to z. 

To conclude this subsection, let us say that matrices 
like eij or -qij have in practice a finite dimension equal to 
Ng X iVc, where Ng is number of reciprocal lattice vec- 
tors considered in the numerical calculations. It is also 
worth stressing that the simple Fourier factorization used 
above for the products like etjEj, which is exact when 
Ng — >■ cxD, may lead in some cases to serious convergence 
problems when truncating the matrices iij. The reason 
is that both the permittivity tensor and the electric field 
can exhibit discontinuities at the interfaces between dif- 
ferent materials. The correct Fourier factorization of this 
type of products when Ng is finite is discussed in detail 
in Appendix A. 



B. Band structure of a single layer 

Now our task is to solve the Maxwell's equations in mo- 
mentum space derived in the previous subsection for the 
case of an unbounded layer. In this case, the fields have 
a z dependence typical from plane waves, i.e. exjp(iqz). 
Moreover, the H field will be expanded in basis set of 
zero divergence, to satisfy V • H = 0, and the coefficients 
in this expansion will be determined by substituting into 
Maxwell's equations. 

Following Ref. [2^, the magnetic field is expanded in 
terms of z propagating plane waves as follows 



H(r,z) 



E 

G 



= (G) 



1 



{kx + Gx)2 



;(G) 



y---{ky + Gy)z 



(15) 

i(k+G)-r+igz 



where x, y, and z are the Cartesian unit vectors and 
(j)xiG) and (f)y{G) are the expansion coefficients to be 
determined. Notice that this expression satisfies V • 
H = 0. Now, it is convenient to rewrite the previ- 
ous expression in momentum representation. By defin- 
ing the vectors (f>x = [4>x{Gi), 4)x{G2), ■ . .]'^ and (j)y = 
[0j,(Gi), (/)j,(G2), • ■ • ]"^, we can write 



Hz) 



bx^ + (t>yy - ^{kx(l)x + ky(l)y)z \ , (16) 



where kx and ky are the diagonal matrices defined in 
Sec. II. For what follows, it is convenient to rewrite this 
last equation in the following vector notation 

h{z) = e''^'' (^(l>x,(l>y,-^{kx(l)x + ky't>y)j ' (17) 

where let us recall that every entry in this column vector 
is a vector of dimension Nq. With this vector notation, 
Eqs. ([9l fTT|) can now be written as 



Ch[z) — ee(z), 
where the block matrices C and e 



C 



q\ -ky 
-qi 6 kx 
ky kx 



^xx ^xy ^xz 
iyx f-yy ^yz 



(18) 



(19) 



^zx ^zy ^zz 

On the other hand, Eqs. (|12m4p adopt now the form 

C'^e{z) ^oj^h{z). (20) 

From Eq. (|18p we obtain the following expression for the 
electric field in momentum representation 



e{z) = fjCh{z), 



(21) 



where 77 = e^^. Substituting this expression in Eq. (j20p 
we obtain the following closed equation for the magnetic 
field in momentum representation 



C^fiCh{z) = ijj^h{z), 



(22) 



which defines an eigenvalue problem for lo^ . Indeed, only 
two of the three identities obtained from this equation, 
one for each x, and z, are independent. From the first 
two identities, and using Eq. (IT7|) . we obtain the following 
equations determining the allowed values for q 



^29^ + Ag + A + ^-1^ ) = 0, (23) 



where = 
defined by 



^x J V^y 



and the 2x2 block matrices An are 



A2 
Ao 



A, 



Vyy 

^Vxy 

(a) 



-A, 



(b) 



Ai = A'( 



(a) 



Af^ 



kyTjz 



kyVzy 


kyTjzX 


kx Vzy 


kx ^ZX 




Vyyk-xk 


"x / 


V Vxxkyk 



Vyzky fjyzkx 

^Xzky ^Xzkx 



/ kyTfzxkykx kyfjzykxkx kyfjzxkyky kyTfzykxky \ 

\ kxfjzykxkx kxfjzxkykx kxTjzykxky kxfjzxkyky J 



1 
1 



(24) 
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In general, Eq. (|23p is a so-called rational eigenvalue 
problem. This problem belongs to the category of non- 
linear eigenvalue problems, which continues to be a chal- 
lenge in the field of numerical analysis. However, we 
have found that a simple linearization strategy allows us 
to solve such an eigenvalue problem in all the examples 
that we have studied. The details of this method are ex- 
plained in Appendix B. It is worth stressing that so far 
the scattering approach has only been applied to situa- 
tions where the materials are isotropic^^ or in cases in 
which the magneto-optical activity is such that the off- 
diagonal components of the permittivity tensor involving 
the z component are zero^^ (txz = £yz =0). In those 
cases, Eq. ((23| reduces to 

= -A2q^(l), (25) 

which is a generalized eigenvalue problem for q^, which 
can be solved with standard techniques of linear algebra. 
Notice that, as explained in the introduction, those cases 
exclude, for instance, the analysis of the Kerr effect in 
the transversal configuration. 

The solution of Eq. (1^51) provides 47Vg non-vanishing 
complex eigenvalues for q. Half of these eigenvalues lie 
in the upper half of the complex plane and half of them 
in the lower half. Finally, let us say that in the case 
of spatially uniform slabs, Eq. (1^^ reduces to a quar- 
tic equation for q, which is well-known in the context of 
wave propagation in anisotropic media. This is shown in 
Appendix C. 

C. Electric and magnetic field 

The next step toward the complete solution of the 
Maxwell's equations in a multilayer structure is the de- 
termination of the fields in the different layers, which can 
be done by expressing them in terms of the propagating 
wave eigenstates defined in the previous subsection. To 
be precise, the fields can be expressed as a combination 
of forward and backward propagating waves with wave 
numbers and complex amplitudes a„ and 6„, respec- 
tively. These amplitudes will be later on determined by 
using the boundary conditions at the interfaces and sur- 
faces of the multilayer structure. Since the boundary 
conditions are simply the continuity of the in-plane field 
components, we shall focus here on the analysis of the 
field components Cx, e^, hx, and hy. 

From the momentum representation of H in Eq. (IT71) , 
the in-plane components of h can be expanded in terms 
of propagating waves as follows 




where d is the thickness of the layer. Here, a„ is the co- 
efhcient of the forward going wave at the z = interface. 



and bn is the backward going wave at z ~ d. On the other 
hand, (/„ correspond to the eigenvalues of Eq. (l23l) with 
Im{g„} > and Pn are the eigenvalues with Im{p„} < 0. 
Notice also that, contrary to the case of isotropic materi- 
als, here the eigenfunctions of the forward and backward 
propagating waves are different in general. 

To make the notation more compact, we now define 
two 2Ng X 2Ng matrices and <&_ whose columns 
are the vectors 0„ and ipn, respectively. Moreover, 
we define the diagonal 2Ng x 2Ng matrices i+{z) and 
f_(d— z), such that [f+(2;)]„„ = e**^"^ and [L(d— z)]„„ — 
^-ip„{d-z) ^ and the 2A'^G-dimensional vectors h\\{z) — 

[hx{z), hy{z)Y' , a = (ai, a2, . . . )^, and 6 = (6i, 62, . . . )^. 
In terms of these quantities, the in-plane magnetic-field 
components become 

/i||(z) = $+f+(z)a + $_L(d-z)6. (27) 

Similarly, using the momentum representation of E 
from Eq. (PT|) it is straightforward to show that the 
in-plane components of the electric field, e||(z) — 
[—ey{z),ex{z)]'^ (note the skew), are given by 

e||(z) = (^A^a'^^+q-'+Af^^++A2^+q)hiz)a 
+ (^'^^-p-i + Af ^<^^ + A2^-p) L(d - z)b, (28) 

where the ^'s are defined in Eq. ([M)) and we have defined 
the 2Ng X 2Ng diagonal matrices q and p such that qnn = 
g„ and p„„ = p„. 

We can now combine Eq. (P7|) and (|28p into a single 
expression as follows 

f Mil M12 A / f+(z)a \ 
\M21 M22 J \ L{d-z)b J ' 

where the 2Ng x 2Ng matrices are defined as 

Mil - + 4^'*+ + ^2$+g, 

M12 = A^o^<P-p-^ +a[^'^<P-+A2^-p, 
M21 = $+, M22^^-. (30) 

D. The scattering matrix 

The final step in our calculation is to use the scatter- 
ing matrix (^-matrix) to compute the field amplitudes 
needed to describe the different relevant physical quan- 
tities. This part of the calculation is practically inde- 
pendent of the type of materials present in the structure 
(isotropic or anisotropic) and it is nicely explained in sec- 
tion V of Ref. [2^. We just include a brief discussion here 
to make this work more self-contained. 

By definition, the ^-matrix relates the vectors of the 
amplitudes of forward and backward going waves, a/ and 
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bi, where / now denotes the layer, in the different layers 
of the structure as follows 



5*11 512 \ / ai' 
S21 S22 J \ h 



(31) 



The field amplitudes in two consecutive layers are re- 
lated via the boundary conditions for the fields, namely 
the continuity of the in-plane components of the fields in 
every interface and surface. If we consider the interface 
between the layer I and the layer l+l, the corresponding 
boundary conditions read 



l|(0) 
ll(0) 



(32) 



i+i 



where di is the thickness of layer /. From this condition, 
together with Eq. (|29l) . it is easy to show that the am- 
plitudes in layers / and / + 1 are related by the interface 
matrix /(/, Z -|- 1) = AIj^^Mi+i in the following way 



fi.+ai 
bi 



I{l,l + l) 



hi I12 
I21 I22 



ai+i 
fi+i-bi+i 

ai+i 
fi+i,-bi+i 



(33) 



where = iL+{di) and /i+i _ = f;+i^_ (dz+i). 

Now, with the help of the interface matrices, the S- 
matrix can be calculated in an iterative way as follows. 
The matrix S{l',l + 1) can be calculated from S{l',l) 
using the definition of S{l',l) in Eq. ([3T|) and the inter- 
face matrix + 1). Eliminating ai and bi we obtain 
the relation between a;/, &;/ and aj+i, from which 

S{1' , I + 1) can be constructed. This reasoning leads to 
the following iterative relations 

Sn{l',l + 1) = [hi- fL+Si2{l',l)l2lY' fL+Sn{l',l) 
Si2il'J+l) = [hi- fL+Sl2{l',l)l2l]~' 

X (/i,+ 5l2(?', 0^22-/12) + 

S2l{l',l + 1) = S22{l',l)hlSll{l\l + l) + S2l{l',l) 
S22{l'.l + l) = S22{l'J)hlSl2il'J + l) + 

822(1', I)l22fl+U-. (34) 

Starting from S{l',l') = 1, one can apply the previous 
recursive relations to a layer at a time to build up S{1' ,1). 

From the knowledge of the S'-matrix one can compute 
all the field amplitudes needed to describe a physical sit- 
uation. Thus for instance, labeling the surface I = and 
the substrate I — N, the calculation of the reflectivity 
and the transmission coefficients requires the knowledge 
of the amplitudes bo and cat, which can be calculated 
from S{0,N). On the other hand, it may be interesting 
to calculate the fields inside the structure, for which we 
need the amplitudes a; and bi . These can be obtained by 



calculating ^(O, I) and S{1, N), and using Eq. (pT|) to get 

ai = [l-Si2{0,l)S2i{l,N)]-' 

X [Sii{0,l)ao + Si2{0J)S22il,N)bN] 

bi = [l-S2i{l,N)Si2{0J)r' 

X [52i(/, ^)5ii(0, l)ao + 822(1, N)bN] . (35) 

III. TMOKE IN PERFORATED IRON FILMS 

In this section we shall illustrate the method just de- 
scribed by analyzing the experiment reported in Ref. [l3 
in which the transverse magneto-optical Kerr effect 
(TMOKE) was studied in a periodically perforated Fe 
film. The TMOKE consists in an intensity change of the 
p-component of the refiected light upon application of a 
magnetic field perpendicular to the plane of incidence of 
the light. In the case of ferromagnetic materials, the 
magnetic field is used to reverse the magnetization M 
of the medium and the TMOKE is characterized by the 
following quantity 



TMOKE : 



Rppi+M) - Rppi-M) 
Rpp{+M)+Rpp{-M)' 



(36) 



where i?pp(±M) are the refiectivity along the p-channel 
for the two opposite magnetizations, which in this con- 
figuration are perpendicular to the incidence plane and 
parallel to the layers of the structure. As explained in 
the introduction, this arrangement induces off-diagonal 
components of the permittivity tensor of the ferromag- 
netic material in the z-direction (direction of the growth 
of the multilayer) and therefore, it requires the use of the 
generalized scattering approach described in the previous 
section. 

The structure studied in Ref. [l3 is described schemat- 
ically in Fig. [1] It consists of a Fe film (100 nm thick), 
which is perforated with a periodic array of subwave- 
length circular holes (diameter of 297 nm) forming a 
triangular lattice with a lattice parameter of 470 nm. 
The Fe film was prepared on a Si(lll) substrate and the 
structure contains additionally a seed layer of Ti (2 nm 
thick) and a capping layer of Au (2 nm thick), which 
were included to form a smooth Fe film and to prevent 
a subsequent oxidation of the surface, respectively. In 
our calculations we used the energy dependent permit- 
tivities taken from ellipsometric measurements of 20 nm- 
thick continuous films, and the off-diagonal elements of 
the ferromagnetic material have been extracted from Po- 
lar Kerr measurements (both rotation and ellipticity) as 
described in Ref. Hf]. Let us emphasize that in the case 
of the Fe film, the permittivity tensor in the transversal 
configuration described in Fig. [TJa) adopts the following 
form'^'^ 



(37) 
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FIG. 1: (Color online) (a) Schematic representation of the 
system under study, which consists of perforated Fe films with 
a periodic array of circular holes forming a triangular lattice. 
Here, one can see the angle definitions and the geometrical 
parameters of the hole array, (b) The triangular lattice both 
in real and in reciprocal space, and definition of the basis 
vectors. 
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FIG. 2: (Color online) Experimental and theoretical results 
for the reflectivity of the demagnetized structure along the p- 
channel, Rpp, as a function of the wavelength of the incident 
light. As indicated in the panels, the results are shown for two 
different high symmetry crystallographic directions ip = 0° 
and tfi = 30° and for various angles of incidence 8. 



where e is the permittivity function of the non- 
magnetized film and Cxz = aM, where M is the mag- 
nitude of the magnetization at saturation. Let us recall 
that in the transverse configuration M — My, i.e. the 
magnetization is parallel to the Fe film and perpendicu- 
lar to the plane of incidence. 

Since the holes of our structure are circular, the Fourier 
expansion of the permittivity, Eq. (|5]), can be calculated 
analytically)^ For holes of radius r and permittivity com- 
ponents in a material with permittivity e^, we have 



,(G) = 



2(e^,-e|)/3^i(Gr)/(Gr) 



+ /3(e 



if G 7^0 
if G = 0, 



(38) 



where /3 is the fraction of the area occupied by the 
holes, and Ji is a Bessel function of first kind. In 
the case of a triangular lattice with lattice constant ao, 
/? = {2lV?>W/al. 

In the numerical calculations performed to obtain the 
results that we are about to describe we have truncated 
the Maxwell equations by setting a high-momentum cut- 
off and we have employed the fast Fourier factorization 
described in Appendix A. In particular, the results shown 
in what follows were obtained by using Nq = 367 lattice 
vectors, which suffices to converge the different physical 
properties discussed here (see Appendix A). 

Let us start our discussion of the results by describ- 
ing the reflectivity in this multilayer system when the 
Fe film is demagnetized. In the upper panels of Fig. [5] 
we reproduce the experimental results for p-polarized 



light obtained for two different high symmetry crystal- 
lographic directions </? = 0° and Lp ~ 30° and for various 
angles of incidence 6 (see Fig. [Ua) for a definition of 
these angles))^ The most prominent feature is the ap- 
pearance of a dip which is red shifted as the angle of 
incidence 9 is increased. Notice that the red shift de- 
pends on the crystallographic direction, and it is more 
pronounced for ip = 0°. Such a feature is absent in the 
case of s-polarized light (not shown here) and it can be 
attributed to the excitation of surface plasmon polaritons 
(SPPs), as we shall discuss below. In the lower panel of 
Fig. [5J we show the corresponding theoretical results cal- 
culated with the scattering approach assuming that exz 
in Eq. (j37p is zero. As one can see, our calculations nicely 
reproduce the experimental trends. The theoretical dips 
appear to be more pronounced than in the experiment, 
which we attribute to the unavoidable inhomogeneities 
in the periodic array of holes in the Fe film. 

The corresponding results (both experimental and the- 
oretical) for the TMOKE, as defined in Eq. are 
displayed in Fig. |31 Notice that the theoretical results, 
in good agreement with the experiment, show that the 
TMOKE can be resonantly enhanced at wavelengths that 
follow closely those in which the dips in the reflectivity 
appear. Notice that at resonance the TMOKE signal in- 
creases by roughly a factor of 2 with respect to value 
at off-resonant wavelengths. It is worth stressing that, 
as illustrated in Fig. HI the signal for the continuous Fe 
film (no perforated) is featureless in the spectral range 
considered here. 
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FIG. 5: (Color online) Theoretical results for the reflectivity 
and TMOKE as a function of the wavelength and angle of 
incidence 9 for tp — 0° and ip = 30°. The dashed lines in the 
different panels correspond to the resonant condition for the 
excitation of SPPs, as described by Eq. (|40|) . 



FIG. 3: (Color online) Experimental and theoretical results 
for the TMOKE as a function of wavelength for tp = 0° and 
if — 30° and for various angles of incidence 9. 



films. Thus, ignoring the thin Au layer, which is prac- 
tically transparent, the complex wave vector of the SPP 
modes is given by^ 



In order to understand the origin of the peaks in the 
TMOKE and the corresponding dips in the reflectivity, 
we have investigated these quantities in a more system- 
atic way. In Fig. [5] we present the results for these two 
quantities as a function of the wavelength and of the an- 
gle of incidence 9. In this figure we can observe again 
the appearance of the dips in the reflectivity, which are 
accompanied by pronounced peaks in the TMOKE. The 
shape of the TMOKE and the dispersion of the peaks 
with 9 suggest that these features originate from the ex- 
citation of the SPPs of this structure. To confirm this 
impression we have to calculate the matching condition 
for the excitation of these surface modes. For this pur- 
pose, we need first to determine the dispersion relation 
of the SPPs and, as an approximation, we shall assume 
that it is given by the dispersion relation for continuous 
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FIG. 4: (Color onhne) Theoretical resuhs for the TMOKE 
for the non-perforated multilayer structure (formed by uni- 
form slabs) as a function of wavelength for various angles of 
incidence 6. 



>(A) 




(39) 



where A is the light wavelength, e is the permittivity of 
the demagnetized Fe film, and we have used the fact that 
the incidence medium is air. Now, the matching condi- 
tion, which is the conservation of the parallel wave vector, 
can be written as 



|Re {/cgpp (Am ,n2 )} 



(40) 



where k|| = (27r/A) sinfiix is the in-plane wave vec- 
tor of the incoming light (in the transverse configura- 
tion described in Fig. [Tfa)) and G„j„2 is a recipro- 
cal lattice vector, which with our choice for the recip- 
rocal lattice basis vectors (see Fig. [TJb)) is given by 
G„i,„2 = (27r/aoV^){[(2n2 - rii) cosi^ -I- ni\/3sin(^]x -I- 
[ni\/3cosip + {ni — 2n2) sin ip]y} . The condition of 
Eq. (HH)) tells us at which (discrete) wavelengths can the 
SPPs be excited for a given angle of incidence. We have 
solved Eq. (|40l) numerically and found that for = 0° 
the only mode that can be excited in the wavelength 
range analyzed here is Ao,-i, while for if = 30° we have 
two possibilities: Ao,-i and A_i__i, which indeed corre- 
spond to the same wavelength. In Fig.[5]we have included 
as dashed lines the relation between the resonant wave- 
length and the angle of incidence 9 for these two cases. 
As one can see, these relations nicely describe the posi- 
tions of both the dips in the reflectivity and the peaks 
in the TMOKE. This strongly suggests that the resonant 
enhancement of the TMOKE is due to the excitation of 
SPPs in our layered structure. 
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IV. CONCLUSIONS 



Appendix A: Fast Fourier factorization 



We have presented in this work a generahzation of the 
scattering-matrix approach to describe the propagation 
of electromagnetic waves in periodically patterned mul- 
tilayer structures containing materials with any kind of 
optical activity and anisotropy. This generalized formal- 
ism enables us to tackle important physical problems that 
have been traditionally out of the scope of this approach. 
Thus for instance, the method can be applied to describe 
all the basic magneto-optical effects in any possible con- 
figuration in magnetic structures. Moreover, the method 
can also be used to study the wave propagation in peri- 
odic structures containing an arbitrary number of bire- 
fringent/dichroic layers. 

We have illustrated the use and capabilities of the 
method by analyzing a recent experiment in which the 
transverse magneto-optical Kerr effect (TMOKE) was in- 
vestigated in a Fe film with a periodic array of subwave- 
length holeSii'^ We have shown, in excellent agreement 
with the experiment, that the TMOKE signal can be res- 
onantly enhanced when the samples are illuminated with 
an appropriate wavelength, and we have attributed this 
phenomenon to the excitation of surface plasmon polari- 
tons. This resonant enhancement of the magneto-optical 
signal is closely related to the phenomenon of extraordi- 
nary optical transmission (EOT))^ which indeed takes 
place in these perforated Fe films. The systematic anal- 
ysis of the interplay between the EOT phenomenon and 
the different magneto-optical effects in perforated mag- 
netic films will be the subject of a forthcoming publica- 
tion. 
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The scattering approach, as formulated in section |TI1 
is known to have important convergence problems when 
metals are involved. These problems are specially pro- 
nounced when the structures contain noble metals and 
the infrared range is investigated. These problems are 
well-known in the theory of gratings^ and it has been 
understood that they originate from the incorrect fac- 
torization of the product of two periodic discontinuous 
functions. Such a product appears, in particular, in the 
constitutive relation D = eE, where D is the displace- 
ment vector. When calculating the Fourier components 
of D in section Hi Al we have used the so-called Laurent's 
rule. This rule states that the Fourier components /i„ 
of the product h{x) of two arbitrary functions f{x) and 
g{x) are given by 



(Al) 



Although this result is correct, as long as the sum extends 
to infinity, it is not always correct when one truncates the 
series, as we do numerically. This was recognized by Li)22 
who established the following rules for factorization: 

1. Let h{x) = f(x)g(x) and cither f[x) or g(x) be 
continuous at some x ^ xq. The other quantity 
may be discontinuous there. Then, Laurent's rule 
applies, i.e. 



[h] = [[.fm. 



(A2) 



Here, [g] denotes a column vector constructed with, 
let us say, Nq Fourier components g„ and by [[/]] 
we denote the Nq x Nq Toeplitz matrix whose 
(n, m) entry is /„ 

2. Let h{x) = f{x)g{x) and both f{x) and g{x) be 
discontinuous at some x — xq, but the product 
f{x)g{x) be continuous there. Then, the so-called 
inverse rule holds, which is given by 



[h] = 



(A3) 



3. Let h{x) = f{x)g{x) and both f{x) and g{x) be 
discontinuous at some x = xq and the product 
f{x)g{x) be discontinuous there as well. Then, the 
product of the two functions in Fourier space can- 
not be formed by either the Laurent's rule or the 
inverse rule. 

Obviously, in our analysis of the Maxwell's equations 
in section III A[ see Eqs. ([ MTTj) . we are violating these 
factorization rules. We are simply using the Laurent's 
rule, although in the interface between different materials 
we may have concurrent discontinuities in both the per- 
mittivity tensor and the electric field, and in some cases 
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the product (the displacement vector) is discontinuous as 
well. Thus, our goal now is to reformulate the Maxwell 
equations in momentum space in order to respect the fac- 
torization rules stated above. For this purpose, we make 
use of the so-called fast Fourier factorization put forward 
by Popov and Neviere in Ref. [SSl 

For the sake of concreteness, let us consider a two- 
dimensional periodic system consisting of an array of cir- 
cular holes or circular pillars, as in the structure of section 
mil Now, let us define a vector with the continuous com- 
ponents of the E and D fields, i.e. G = [Et, Dn, E^]'^ . 
Here, Et is the tangential component of the electric field 
in the xy plane, £)„ is the normal component of the dis- 
placement vector in the xy plane, and Ez is the z com- 
ponent of the electric field. These three components are 
continuous in the xy plane when we cross the bound- 
ary of a hole (or pillar) and the permittivity tensor un- 
dergoes a discontinuity. Now, let us establish the rela- 
tion between these field components and the three Carte- 
sian components of the electric field G — EE, where 
E = [Ex, Ey, Ez]"^ . There are many possible choices for 
F. We choose to express its matrix elements in terms of 
the polar angle 0(a;, y) defined as re**'^'^' ~ x + iy. It is 
straightforward to show that 

-s c 

1 

where c and s are abbreviations for cos and sin cf), re- 
spectively. 

We now define the inverse of this matrix C = F^^. 
Thus, E = CG. Let us recall the constitutive relation 
D = eE, where e is given by Eq. ([1]). This relation can 
be now written as 



r> = eC ■ G ^ eC ■ EE, 



(A5) 



whose elements are now expressed as the product of a 
discontinuous function and a continuous one. Thus, using 



Laurent's rule for the first term of the product and the 
inverse rule for the second one, the Fourier components 
of the displacement vector can be calculated as 



(A6) 



This indicates that the Toeplitz matrix of the index ten- 
sor 77 in the formalism of section has to be calculated 
as follows 



(A7) 



This is indeed the only change that we need to introduce 
in the formalism to improve significantly the convergence 
in the problematic cases. Notice that in practice this 
requires the calculation of the Toeplitz matrix of several 
trigonometric functions, which in general has to be done 
numerically. 

For the sake of completeness, we now provide simplified 
expressions for [[77]] in some cases of special interest for us. 
First, in the case of an isotropic material, for which e = 
el, it is straightforward to show that Eq. (jA7p reduces 




(A8) 



-' + [[X]][[c']] [[X]][[cs]] 

miles]] [[l/e]]-[[X]][[c^]] 

M' 



where [[X]] — [[1/e]] — [[e]]""'". This requires, in partic- 
ular, the calculation of the Fourier components of the 
trigonometric functions cos^ (j) and cos (j) sin 0, which can 
be easily done numerically. Notice that these compo- 
nents are independent of the wavelength of the light and 
therefore, they can be calculated once and forever for a 
given structure. On the other hand, for the description 
of the TMOKE we have to consider a permittivity tensor 
given by Eq. (|37p . In this case, the Fourier components 
of the index tensor are given by 



J 



^]] + m][[c']] [[Y]][[cs]] 

m][[cs]] [[l/e]]-'-[[Y]][[c^]] 
ex.]]-[[Z]][[c']] -[[Z]][[cs]] 



[exz]] + [[Z']][[c']- 
[[Z']] [[cs]] 

M] + [[W]]M] 



(A9) 



r 



where 

[[Y]] = [[1/6]]-' -[H], 

[[Z]] - [[e.^/e]] [[1/e]]-' [[exz]], 

[[Z']] = [[l/e]]-'[[exje]]-[[exz]], 

[[W]] = [[elJe]]~[[exz/e]][[l/e]]-%xJe]].iA10) 

Again, one just needs the evaluation of the Fourier com- 
ponents of both cos^ (j) and cos (j) sin 0. 



It is worth stressing that in the choice of the normal 
vectors entering in the Fourier factorization there is a 
freedom that one can use to further improve the conver- 
gence of the calculations. For a discussion of this issue, 
see Ref. M 

To conclude, we now want to illustrate the conver- 
gence of the results using the fast Fourier factorization 
described in this appendix. In Fig. |5] we show an ex- 
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wavelength (nm) wavelength (nm) 



We have thus converted the problem into a generalized 
linear eigenvalue problem that can be solved with stan- 
dard linear algebra techniques. The obvious disadvantage 
of this simple procedure is that one increases the dimen- 
sion of the problem by a factor of 3. In this sense, it 
may be advantageous in some cases to implement other 
methods like, for instance, the iterative Newton method. 
In any case, and as illustrated in the previous appendix, 
we have not found problems to converge the calculations 
detailed in this work with the linearization procedure. 



FIG. 6: (Color online) Reflectivity in the absence of magnetic 
field (a) and TMOKE (b) for the multilayer structure studied 
in section IIIII for 9 = 25° and <^ = 0° as a function of the 
wavelength of the incident light. The difi'erent curves corre- 
spond to results obtained with different number of reciprocal 
lattice vectors No- 



ample of the results obtained for the reflectivity and the 
TMOKE for the structure studied in section Hill In this 
figure, the different curves correspond to different val- 
ues of Nq^ which is the number of reciprocal lattice vec- 
tors taken into account in the calculations upon setting 
a high-momentum cutoff. As one can see, it is possible 
to converge the calculations to a high precision in the 
whole range of wavelengths. Moreover, the convergence 
is rapid and uniform. It is important to emphasize that 
in order to get results of similar quality for this example 
without the use of the fast Fourier factorization, values of 
Nq even larger than 1000 are required (not shown here). 

Appendix B: Solving the nonlinear eigenvalue 
problem 

We detail here a simple strategy to solve numerically 
the nonlinear eigenvalue problem of Eq. (j23l) . which we 
have found to work without any problem in all the cases 
that we have considered. The first step is to multiple 
both sides of Eq. by q to convert it into the following 
cubic eigenvalue problem 

(S3«' + S2g'+Si(? + 6o)</' = 0, (Bl) 

where we have defined Bn = ^n-i- Now, the simplest 
strategy to solve this cubic problem is to use a standard 
linearization procedure.'*''^ The idea goes as follows. We 
first define the following vectors 

A„=g"-V; n= 1,2,3. (B2) 

From this definition, and using Eq. (jBip . it is easy to 
show that the vectors A„ satisfy the following equation 

1 0\/Ai\ /lO 0\/Ai\ 
1 As ^ g 1 Aa . 

Bo Bl 5-2 J \X3 J \0 -Bs J \X3 J 

(B3) 



Appendix C: Spatially uniform slabs 

A multilayer structure may contain some uniform 
(non-structured) layers. In particular, this is always the 
case for the medium of incidence and for the substrate 
layer. In this sense, it is interesting to discuss how the 
formalism discussed in section |lT] is simplified in the case 
of uniform slabs. In this case, the permittivity tensor is 
diagonal in momentum space: {iij)G,G' = eii(0)(5cj.G'- 
This implies that all the matrices in momentum repre- 
sentation are also diagonal. The eigenvalue problem of 
Eq. (|23l) leads to the following quartic secular equation 
for q{G). Focusing on G = 0, this equation reads 

4 

Y,Dnq" = 0, (CI) 

n=0 

where the coefficients are given by 

-^3 {^xy'Hyz ~t~ VyxVzy ^yy{j]xz ~t~ ^za;)] 

~\~ky [l^yxVxz ~t- ^xyVzx ^xxi'Hyz Vzy)] : 

D2 = k,j. \Tjyy{Tjxx + Vzz) 'Hxy'Hyx VyzVzy] 
+ky [VxxiVyy + Vzz) - 

[ilxz{riyz +11zy) + Vyzillzx - Vxz) 
- VzziVxy + Vyx)] - i^^iVxx +'nyy), 

Di = kl[rixyriy^ + rjy^ri^y - riyy{rix^ + rj^.^)] 

'^^'y [VyxVxz + VxyVzx ~ VxxiVyz + Vzy)] 
-\-kj,ky \j]xy^zx ~t- TjxzVyx ^xx{j]yz ~^ ^^y)] 
-\-kykx [l]yxVzy ~t~ VyzVxy - VyyiVxz + Vzx)] 

+^'^ [kliVxz + Vzx) + kl{riy^ + ri,^y)\ , 

Do = kxiVyyVzz - VyzTlzy) + kyiVxxVzz - Vxzrizx) 

+kx ky 

[VxzVzy + VyzVzx - VzziVxy + Vyx)] 
~^kykx [VyzVzx ~l" TJxzVzy Vzz iVyx+Vxy)] 
-\-kxky [tIzz{i1xx + Vyy) ^" VxyVyx ~ VxzVzx ~ VyzVzy] 
[w^ - kliriyy + 77^2) - ky{T]xx + Vzz) 
-{- kxkyijjxy "I" Vyx )] ■ (C2) 

For G 7^ 0, one just needs to replace kx^y by kx^y + Gx,y 
Eq. (jCip has been previously derived (in terms of the 
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components of the permittivity tensor) in the context of 
the analysis of uniform multilayer structures containing 
magneto-optical and anisotropic materials.— This equa- 
tion simplifies in several limiting cases. Thus for in- 
stance, if we consider the typical configuration for mea- 
suring the TMOKE, then the permittivity tensor is given 
by Eq. p7|) . In this case, D3 — Di ~ and setting 



the solutions of Eq. (jCl|) are 



and ^2 



matrix defined in 
form (for G = 0) 



Eq. 



hi 



Moreover, in this case the layer 
29]) adopts the following simple 



On the other hand, for an isotropic layer e = 
this case D^, = Di = Q and = ew^ — (fc^ - 
corresponding layer matrix reads now (for G 



M = 



i 



-Mil 

i 



el, and in 
-kl). The 
= 0) 



(C4) 



where 1 is the 2x2 unit matrix and 



M 



Vxxq2 ~ Vxz 





V 



kx 





1 







Tlxxq2 Vxzkx 





1 



(C3) 



Mil = 



Tjkx ky 



(C5) 
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